This notebook draws from previous notebooks on relative cumulative frequency analysis and PCA analysis of metabolic cage data.

RMR data before addition of the wheels was analyzed and presented previously.

Dataset

These datasets are from RMR data collected for mice after introduction to the metabolic cages (days 0-6, see image).

overview

A total of 24 mice were run over the length of 2 experiments.

Experiment 1: The first set of mice, run in February

Start of Met Cage run	Metabolic Cage #	Ear Tag	DOB	Gender	Genotype	Body Mass on BMR day (g)
2/4/2020	1	1049	10/08/19	M	WT	27.4
2/4/2020	2	1050	10/08/19	M	L-Bmal1-KO	26.6
2/4/2020	3	1053	10/07/19	M	WT	31.5
2/4/2020	4	1054	10/07/19	M	WT	30.5
2/4/2020	5	1055	10/07/19	M	L-Bmal1-KO	27.5
2/4/2020	6	1109	10/11/19	M	WT	28.3
2/4/2020	9	1110	10/11/19	M	WT	27.9
2/4/2020	10	1111	10/11/19	M	L-Bmal1-KO	26.7
2/4/2020	11	1112	10/11/19	M	WT	28.5
2/4/2020	12	1214	10/21/19	M	L-Bmal1-KO	25.3
2/4/2020	13	1215	10/21/19	M	WT	27.9
2/4/2020	14	1216	10/21/19	M	L-Bmal1-KO	27.6

Experiment 2: The second set of mice, run in March

Start of Met Cage run	Metabolic Cage #	Ear Tag	DOB	Gender	Genotype	Body Mass on BMR day (g)
3/10/2020	1	1555	11/18/19	M	WT	26.4
3/10/2020	2	1692	12/06/19	M	WT	25.6
3/10/2020	3	1693	12/06/19	M	WT	24.7
3/10/2020	4	1695	12/06/19	M	WT	25.6
3/10/2020	5	1699	12/06/19	M	WT	28.7
3/10/2020	6	1556	11/18/19	M	L-Bmal1-KO	26.2
3/10/2020	9	1557	11/18/19	M	L-Bmal1-KO	27.1
3/10/2020	10	1558	11/18/19	M	L-Bmal1-KO	26.9
3/10/2020	11	1559	11/18/19	M	L-Bmal1-KO	24.7
3/10/2020	12	1691	12/06/19	M	L-Bmal1-KO	28
3/10/2020	13	1694	12/06/19	M	L-Bmal1-KO	26.4
3/10/2020	14	1700	12/06/19	M	L-Bmal1-KO	27.7

Create data

Data generated using a python script, called via shell command in R. A conda environment is required (see PCA notes ).

Bmal.py

#!/usr/bin/env python

"""
A python script that uses the mousetrapper library to prepare WT vs KO mice data for import to R.
Must be run from same directory.
"""

import numpy as np
import pandas as pd

import lib.mousetrapper as mt

__author__ = "Sumeed Yoyo Manzoor"
__email__ = "smanzoor@uchicago.edu"
__version__ = "0.0.1"
__experiment__ = "12 WT vs 12 L-Bmal1-KO mice"

datasets = "../datasets/"
data = "../data/"
outfolder = data + "R/"

experiment1_RMR = mt.Experiment.experiment(
    info="""First run with revised protocol, WT vs L-Bmal1-KO""",
    RT_data=datasets + "Bmal_WT_vs_KO_SPF_RMR_Wheel_2020-02-10_rt.csv",
    macro13_data=datasets + "bmal_wt_vs_ko_spf_rmr_wheel_fixedwheel_2020-02-10_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv", 
    t_ctrl={"SPF WT":[1,3,4,6,9,11,13]},
    t_exp={"SPF L-Bmal1-KO":[2,5,10,12,14]},
    groups={
        "SPF WT":[1,3,4,6,9,11,13],
        "SPF L-Bmal1-KO":[2,5,10,12,14]
    },
    meta=["WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","L-Bmal1-KO","WT","L-Bmal1-KO"],
    covar=[27.4,26.6,31.5,30.5,27.5,28.3,27.9,26.7,28.5,25.3,27.9,27.6]
    )

experiment2_RMR = mt.Experiment.experiment(
    info="""Second run with revised protocol, WT vs L-Bmal1-KO""",
    RT_data=datasets + "Bmal_SPF_KF_RMR_Wheel_2020-03-16_16_31_rt.csv",
    macro13_data=datasets + "bmal_spf_kf_rmr_wheel_2020-03-16_16_31_fixedwheel_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv",
    t_ctrl={"SPF WT":np.arange(1,6).tolist()},
    t_exp={"SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()},
    groups={
        "SPF WT":np.arange(1,6).tolist(),
        "SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()
    },
    meta=np.array(["WT","L-Bmal1-KO"]).repeat([5,7]).tolist(),
    covar=[26.4,25.6,24.7,25.6,28.7,26.2,27.1,26.9,24.7,28,26.4,27.7]
    )

experiment1_RMR.macro13_lightcycles(days=2)
experiment2_RMR.macro13_lightcycles(days=2)

experiment1_RMR.macro13_trimdays(days=2, timeformat=False)
experiment2_RMR.macro13_trimdays(days=2, timeformat=False)

experiment1_RMR.macro13_dark.to_csv(outfolder + "1_dark.csv", index=False)
experiment1_RMR.macro13_light.to_csv(outfolder + "1_light.csv", index=False)
experiment1_RMR.macro13_trimmeddays.to_csv(outfolder + "1_total.csv", index=False)
pd.DataFrame({"covar":experiment1_RMR.covar}).to_csv(outfolder + "1_covar.csv", index=False)

experiment2_RMR.macro13_dark.to_csv(outfolder + "2_dark.csv", index=False)
experiment2_RMR.macro13_light.to_csv(outfolder + "2_light.csv", index=False)
experiment2_RMR.macro13_trimmeddays.to_csv(outfolder + "2_total.csv", index=False)
pd.DataFrame({"covar":experiment2_RMR.covar}).to_csv(outfolder + "2_covar.csv", index=False)

shell("conda activate met_cages_development && cd py && python Bmal_wheel.py")

Import data

reqpkg <- c("magrittr","dplyr","ggplot2","ggpubr","DT","reshape2")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "DT"
## [1] '0.13'
## [1] "reshape2"
## [1] '1.4.3'

theme_set(theme_pubr())
select <- dplyr::select

path <- "data/R/"

exp1 <- read.csv2(paste0(path, "1_total.csv"), sep = ",") %>% 
  set_colnames(paste0("exp1_",colnames(.)))
exp1_dark <- read.csv2(paste0(path, "1_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("exp1_",colnames(.)))
exp1_light <- read.csv2(paste0(path, "1_light.csv"), sep = ",") %>% 
  set_colnames(paste0("exp1_",colnames(.)))
exp1_covar <- read.csv2(paste0(path, "1_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

exp2 <- read.csv2(paste0(path, "2_total.csv"), sep = ",") %>% 
  set_colnames(paste0("exp2_",colnames(.)))
exp2_dark <- read.csv2(paste0(path, "2_dark.csv"), sep = ",") %>%
  set_colnames(paste0("exp2_",colnames(.)))
exp2_light <- read.csv2(paste0(path, "2_light.csv"), sep = ",") %>% 
  set_colnames(paste0("exp2_",colnames(.)))
exp2_covar <- read.csv2(paste0(path, "2_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

covar <- c(exp1_covar[c(1,3,4,6,7,9,11),], exp2_covar[1:5,], exp1_covar[c(2,5,8,10,12),], exp2_covar[6:12,])

total2 <- exp1 %>% 
  bind_cols(exp2) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

light <- exp1_light %>% 
  bind_cols(exp2_light) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

dark <- exp1_dark %>% 
  bind_cols(exp2_dark) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

total <- bind_rows(light, dark)

Datatables

Light

Hide

Select a tab above to view

VO2

light %>% select(contains("VO2")) %>% datatable(filter="bottom")

RQ

Note some data has errors (greater than 1 or below 0.5). These were likely measurement errors.

light %>% select(contains("RQ")) %>% datatable(filter="bottom")

BodyMass

light %>% select(contains("BodyMass")) %>% datatable(filter="bottom")

Food

light %>% select(contains("Food")) %>% datatable(filter="bottom")

Water

light %>% select(contains("Water")) %>% datatable(filter="bottom")

PedMeters

light %>% select(contains("PedMeters")) %>% datatable(filter="bottom")

WheelMeters

light %>% select(contains("WheelMeters")) %>% datatable(filter="bottom")

Dark

Hide

Select a tab above to view

VO2

dark %>% select(contains("VO2")) %>% datatable(filter="bottom")

RQ

dark %>% select(contains("RQ")) %>% datatable(filter="bottom")

BodyMass

dark %>% select(contains("BodyMass")) %>% datatable(filter="bottom")

Food

dark %>% select(contains("Food")) %>% datatable(filter="bottom")

Water

dark %>% select(contains("Water")) %>% datatable(filter="bottom")

PedMeters

dark %>% select(contains("PedMeters")) %>% datatable(filter="bottom")

WheelMeters

dark %>% select(contains("WheelMeters")) %>% datatable(filter="bottom")

Graphs

select_WT <- function(df) {
  d1 <- df %>% select(contains("exp1"))
  for (name in colnames(d1)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(2,5,10,12,14)) {
      d1 %<>% select(-all_of(name))
    }
  }
  
  d2 <- df %>% select(contains("exp2"))
  for (name in colnames(d2)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) > 5) {
      d2 %<>% select(-all_of(name))
    }
  }
  
  d <- bind_cols(d1, d2)
  return(d)
}

select_KO <- function(df) {
  d1 <- df %>% select(contains("exp1"))
  for (name in colnames(d1)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,3,4,6,9,11,13)) {
      d1 %<>% select(-all_of(name))
    }
  }
  
  d2 <- df %>% select(contains("exp2"))
  for (name in colnames(d2)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) <= 5) {
      d2 %<>% select(-all_of(name))
    }
  }
  
  d <- bind_cols(d1, d2)
  return(d)
}

graph_groups <- function(data, channel, title="", type="m") {
  df <- data %>% select(contains(channel))
  df.observation <- 1:nrow(df)
  c <- ncol(df)
  WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
  KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
  df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
  df$observation <- rep(df.observation, c)
  
  if (type=="m") {
    df %>% group_by(variable, observation) %>% 
      summarise(value=mean(value)) %>% 
      ggplot(aes(x=observation, y=value, color=variable)) +
        geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
        geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
        geom_line(size=1) +
        scale_color_manual(values = c("#103EFB", "#FC2A1C")) +
        ggtitle(title) +
        ylab(channel)
  } else if (type=="s") {
    ggplot(df, aes(x=observation, y=value, color=variable)) + 
      ggtitle(title) + 
      geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      scale_color_manual(values = c("#103EFB", "#FC2A1C")) +
      stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", alpha = 0.7) +
      ylab(channel)
  }
}

graph_mice <- function(data, channel, title="") {
  df <- data %>% select(contains(channel))
  df.observation <- 1:nrow(df)
  c <- ncol(df)
  WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt()
  KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt()
  df <- bind_rows(WT, KO) %>% set_colnames(c("mouse", "channel"))
  df$observation <- rep(df.observation, c)
  ggplot(df, aes(x=observation, y=channel, color=mouse)) + 
    geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
    geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
    geom_line() +
    ylab(channel)
}

There are 3 graphs per variable, a graph with the average of the mice in each group, a graph with the confidence limits for the population means (basically standard error), and a graph with each individual mouse plotted.

Hide

Select a tab above to view

VO2

graph_groups(total2, "VO2", title="VO2 by group", type = "m")

graph_groups(total2, "VO2", title="VO2 by group", type = "s")

graph_mice(total2, "VO2", title="VO2 by mouse")

RQ

total2 %>% filter_at(vars(contains("RQ")), all_vars(. < 1 & . > 0.5)) %>% graph_groups("RQ", title="RQ by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

total2 %>% filter_at(vars(contains("RQ")), all_vars(. < 1 & . > 0.5)) %>% graph_groups("RQ", title="RQ by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

total2 %>% filter_at(vars(contains("RQ")), all_vars(. < 1 & . > 0.5)) %>% graph_mice("RQ", title="RQ by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

BodyMass

graph_groups(total2, "BodyMass", title="BodyMass by group", type = "m")

graph_groups(total2, "BodyMass", title="BodyMass by group", type = "s")

graph_mice(total2, "BodyMass", title="BodyMass by mouse")

Food

graph_groups(total2, "Food", title="Food by group", type = "m")

graph_groups(total2, "Food", title="Food by group", type = "s")

graph_mice(total2, "Food", title="Food by mouse")

Water

graph_groups(total2, "Water", title="Water by group", type = "m")

graph_groups(total2, "Water", title="Water by group", type = "s")

graph_mice(total2, "Water", title="Water by mouse")

PedMeters

graph_groups(total2, "PedMeters", title="PedMeters by group", type = "m")

graph_groups(total2, "PedMeters", title="PedMeters by group", type = "s")

graph_mice(total2, "PedMeters", title="PedMeters by mouse")

PedSpeed

graph_groups(total2, "PedSpeed", title="PedSpeed by group", type = "m")

graph_groups(total2, "PedSpeed", title="PedSpeed by group", type = "s")

graph_mice(total2, "PedSpeed", title="PedSpeed by mouse")

AllMeters (R)

graph_groups(total2, "AllMeters_R", title="AllMeters_R by group", type = "m")

graph_groups(total2, "AllMeters_R", title="AllMeters_R by group", type = "s")

graph_mice(total2, "AllMeters_R", title="AllMeters_R by mouse")

AllMeters (M)

graph_groups(total2, "AllMeters_M", title="AllMeters_M by group", type = "m")

graph_groups(total2, "AllMeters_M", title="AllMeters_M by group", type = "s")

graph_mice(total2, "AllMeters_M", title="AllMeters_M by mouse")

WheelMeters

graph_groups(total2, "WheelMeters", title="WheelMeters by group", type = "m")

graph_groups(total2, "WheelMeters", title="WheelMeters by group", type = "s")

graph_mice(total2, "WheelMeters", title="WheelMeters by mouse")

WheelSpeed

graph_groups(total2, "WheelSpeed", title="WheelSpeed by group", type = "m")

graph_groups(total2, "WheelSpeed", title="WheelSpeed by group", type = "s")

graph_mice(total2, "WheelSpeed", title="WheelSpeed by mouse")

ZBreaks

graph_groups(total2, "ZBreak", title="ZBreak by group", type = "m")

graph_groups(total2, "ZBreak", title="ZBreak by group", type = "s")

graph_mice(total2, "ZBreak", title="ZBreak by mouse")

H20

graph_groups(total2, "H2O", title="H2O by group", type = "m")

graph_groups(total2, "H2O", title="H2O by group", type = "s")

graph_mice(total2, "H2O", title="H2O by mouse")

VOC

graph_groups(total2, "VOC", title="VOC by group", type = "m")

graph_groups(total2, "VOC", title="VOC by group", type = "s")

graph_mice(total2, "VOC", title="VOC by mouse")

H2

graph_groups(total2, "H2_", title="H2 by group", type = "m")

graph_groups(total2, "H2_", title="H2 by group", type = "s")

graph_mice(total2, "H2_", title="H2 by mouse")

RCF

The idea behind this was to apply a method to study large metabolic datasets quantitatively. Based on the Gordon paper, we can use relative cumulative frequency to find frequency of data binned by VO2 or RQ reading levels, then calculate EC₅₀. We can use EC₅₀ values for every mouse in a simple t-test or ANCOVA to compare groups. In other words, this can be thought of as a kind of weighted average.

Halatchev IG, O’Donnell D, Hibberd MC, Gordon JI. Applying indirect open-circuit calorimetry to study energy expenditure in gnotobiotic mice harboring different human gut microbial communities. Microbiome. 2019;7(1):158. Published 2019 Dec 12. doi:10.1186/s40168-019-0769-4
Riachi M, Himms-Hagen J, Harper ME. Percent relative cumulative frequency analysis in indirect calorimetry: application to studies of transgenic mice. Can J Physiol Pharmacol. 2004;82(12):1075‐1083. doi:10.1139/y04-117

Relative cumulative frequency

This is a method where data is first ordered by increasing value, then, for every increase along a set increment, a frequency is recorded. Finally, frequency is cumulated. This is all accomplished by the function rcf. In order to run this, however, the range of data needs to be checked to determine suitable range and increments. find_range can be run to check the range on a vector, and will return a recommended range and increment value to submit to rcf.

EC₅₀

By fitting frequencies to intervals, a sigmoidal regression is found. This is dependent on a normal distribution — that is, the mouse should have less frequencies towards the extremes, and the most expression in the middle. For RQ and VO₂, this is true, the mouse will not spend as much time on the extremes of the values. For other data types, the nature and distribution of the data should be considered before using in this analysis.

t-test and ANCOVA

For the Welch t-test, I run an F-test first (var.test), which should be not significant.
ANCOVA is run using the EC₅₀ values as dependent variable and bodyweight (see datasets) as covariate. The ANCOVA function returns a slew of information related to the test. A linear regression is first fit to the dependent variable (DV) and covariate. This may be important (see the Speakman paper ) since it can help elucidate the validity of ANCOVA related to other models, but unfortunately without lean mass data the conclusions to be drawn are limited. Evidently, rarely does energy expenditure regress linearly with total mass. Next, ANOVA is run between DV and covariate. Hopefully, the covariate is not a significant predictor of DV. Then, a boxplot for DV vs independent variable (IV), in this case WT vs KO, and covariate vs IV are shown to check for extreme outliers and IV relationship. Then, basic stats on all variables are displayed, followed by Levene’s test, which should not have a significant result. Sanity checks are run using ANOVA of DV, ANOVA of covariate, then ANOVA with DV normalized to covariate (which should not be used). ANCOVA is finally run. The check_ANCOVA function will pull any significant (<0.05) p values from ANCOVA.
In the future, more complicated models, such as ANCOVA/linear models with more covariates, might be worth implementing.

Libraries

reqpkg <- c("docstring","tibble","scales","purrr","nplr", "car","compute.es","effects","multcomp","pastecs","WRS2")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "docstring"
## [1] '1.0.0'
## [1] "tibble"
## [1] '2.1.3'
## [1] "scales"
## [1] '1.1.0'
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'

Functions

find_range <- function(data, channel, segments=20) {
  data <- data %>% select(contains(channel)) %>% melt()
  r <- range(data$value)
  cut <- max(data$value)-min(data$value)
  out <- list("range"=r, "cut"=cut/segments)
  return(out)
}

rcf <- function(data, channel, id, min, max, b=0.1) {
  #' Relative cumulative frequency function
  #'
  #' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
  #'
  #' @param data numeric vector
  #' @param channel character. Channel to label the rcf output.
  #' @param b double. Number to cut vector.
  #' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
  
  breaks <- seq(min, max, by=b)
  duration.cut <- cut(data, breaks, right = F)
  duration.freq <- table(duration.cut)
  duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
  duration.cumrelfreq = duration.cumfreq/length(data)
  out <- duration.cumrelfreq %>% as.data.frame()
  out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
  out$id <- rep(id, nrow(out))
  return(out)
}

EC50 <- function(x) {
  out <- lapply(x, function(y) {
    nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
  }) %>% 
    sapply(function(z) {
      getEstimates(z, 0.5)$x
    })
  out %<>% set_rownames(NULL)
  return(out)
}

lm_eqn <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                   list(a = format(unname(coef(m)[1]), digits = 2),
                        b = format(unname(coef(m)[2]), digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));
}

lm_eqn_str <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
  return(eq)
}

remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

ANCOVA <- function(df) {
  # Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
  
  cat("linear fit of dependent variable with covariates\n")
  p1 <- ggplot(df, aes_string(x=names(df)[2], y=names(df)[3])) +
    geom_smooth(method = "lm", se = FALSE) +
    geom_point() +
    # geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    ggtitle("Linear model - DV and Covar") +
    xlab("DV") + ylab("Covar")
  print(p1)
  
  m <- lm(df[[3]] ~ df[[2]], df)
  
  summary(m) %>% print()
  lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
  
  plot1 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[2])) +
    geom_boxplot() +
    # ggtitle("DV vs IV") +
    xlab("IV") + ylab("DV")
  plot2 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[3])) +
    geom_boxplot() +
    # ggtitle("Covar vs IV") +
    xlab("IV") + ylab("Covar")
  p2 <- ggarrange(plot1, plot2,
                  labels = c("DV vs IV", "Covar vs IV"),
                  ncol = 2)
  # grid.arrange(plot1, plot2, ncol=2)
  print(p2)
  
  cat("\n_____________________\n")
  cat("Some stats on the variables\n")
  cat("DV - IV:\n")
  by(df[[2]], df[[1]], stat.desc) %>% print()
  cat("Covar - IV:\n")
  by(df[[3]], df[[1]], stat.desc) %>% print()
  cat("levene's test\n")
  leveneTest(df[[2]], df[[1]], center = median) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of DV with IV\n")
  model <- aov(df[[2]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of covariate with IV, independent of DV\n")  
  model <- aov(df[[3]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of DV with IV, normalized to covariate\n")
  c <- df[[2]]/df[[3]]
  model <- aov(c ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANCOVA\n")
  model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
  summary(model) %>% print()
  out <- Anova(model, type= "III")
  print(out)
  return(out)
}

check_ANCOVA <- function(res) {
  if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
  } else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
    print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
  } else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
  }
}

VO2

Hide

Select a tab above to view

Total cycle

total.VO2 <- total %>% select(contains("VO2"))
total.VO2.WT <- total.VO2 %>% select_WT()
total.VO2.KO <- total.VO2 %>% select_KO()

rcf.WT <- lapply(total.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.001))
rcf.KO <- lapply(total.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.001))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})

# overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in complete cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 1.71539878847918"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 1.92269230171977"

Run F & t test

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.5506, num df = 11, denom df = 11, p-value = 0.4787
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4463848 5.3863404
## sample estimates:
## ratio of variances 
##           1.550607

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -1.7594, df = 21.02, p-value = 0.09306
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3565851  0.0297364
## sample estimates:
## mean of x mean of y 
##  1.776285  1.939710

Run ANCOVA

res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()

## linear fit of dependent variable with covariates

## `geom_smooth()` using formula 'y ~ x'

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6789 -0.9095 -0.1433  0.7997  4.0970 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.5857     2.7174   10.52 4.76e-10 ***
## df[[2]]      -0.7256     1.4512   -0.50    0.622    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.654 on 22 degrees of freedom
## Multiple R-squared:  0.01124,    Adjusted R-squared:  -0.03371 
## F-statistic:  0.25 on 1 and 22 DF,  p-value: 0.622
## 
## y = 29 + -0.73x, r^2 = 0.0112

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.58030227   2.30114728   0.72084501 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  23.27651438   1.91601402   1.93970953   0.05816211   0.12801395   0.04059398 
##      std.dev     coef.var 
##   0.20147947   0.10387095 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.36140182   2.11187515   0.75047334 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  21.31542198   1.67948358   1.77628517   0.07242541   0.15940726   0.06294529 
##      std.dev     coef.var 
##   0.25088899   0.14124365 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   0.792 0.3831
##       22               
## _____________________
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.1602 0.16025   3.095 0.0924 .
## Residuals   22 1.1389 0.05177                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0004118 0.0004118   5.162 0.0332 *
## Residuals   22 0.0017550 0.0000798                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.1602 0.16025   2.955  0.100
## df[[3]]      1 0.0001 0.00007   0.001  0.971
## Residuals   21 1.1389 0.05423               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.27674  1  5.1030 0.03464 *
## df[[1]]     0.14572  1  2.6869 0.11607  
## df[[3]]     0.00007  1  0.0013 0.97144  
## Residuals   1.13886 21                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

Light cycle

light.VO2 <- light %>% select(contains("VO2"))
light.VO2.WT <- light.VO2 %>% select_WT()
light.VO2.KO <- light.VO2 %>% select_KO()

rcf.WT <- lapply(light.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.001))
rcf.KO <- lapply(light.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.001))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in light cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 1.38854018713538"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 1.44765116831737"

Run F & t test

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 0.44565, num df = 11, denom df = 11, p-value = 0.1959
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1282914 1.5480390
## sample estimates:
## ratio of variances 
##          0.4456457

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -1.445, df = 19.18, p-value = 0.1646
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.24328399  0.04448241
## sample estimates:
## mean of x mean of y 
##  1.373416  1.472816

Run ANCOVA

res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()

## linear fit of dependent variable with covariates

## `geom_smooth()` using formula 'y ~ x'

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4161 -1.0136 -0.0008  0.6401  4.0365 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.484      2.696   8.340 2.94e-08 ***
## df[[2]]        3.341      1.881   1.776   0.0896 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.556 on 22 degrees of freedom
## Multiple R-squared:  0.1254, Adjusted R-squared:  0.0856 
## F-statistic: 3.153 on 1 and 22 DF,  p-value: 0.08962
## 
## y = 22 + 3.3x, r^2 = 0.125

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.21686523   1.92189357   0.70502834 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  17.67379688   1.43491951   1.47281641   0.05721110   0.12592078   0.03927732 
##      std.dev     coef.var 
##   0.19818505   0.13456195 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.17590660   1.51816084   0.34225424 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  16.48098736   1.39123252   1.37341561   0.03819224   0.08406055   0.01750377 
##      std.dev     coef.var 
##   0.13230180   0.09633049 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2399 0.6291
##       22               
## _____________________
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0593 0.05928   2.088  0.163
## Residuals   22 0.6246 0.02839               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0001871 1.871e-04    6.57 0.0177 *
## Residuals   22 0.0006264 2.847e-05                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0593 0.05928   2.644 0.1189  
## df[[3]]      1 0.1537 0.15366   6.852 0.0161 *
## Residuals   21 0.4709 0.02243                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.00023  1  0.0101 0.92079  
## df[[1]]     0.12721  1  5.6728 0.02677 *
## df[[3]]     0.15366  1  6.8522 0.01608 *
## Residuals   0.47093 21                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

## [1] "Varies significantly with DV and Covar,"

Dark cycle

dark.VO2 <- dark %>% select(contains("VO2"))
dark.VO2.WT <- dark.VO2 %>% select_WT()
dark.VO2.KO <- dark.VO2 %>% select_KO()

rcf.WT <- lapply(dark.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.001))
rcf.KO <- lapply(dark.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.001))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in dark cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 2.24612143206034"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 2.59513184736794"

Run F & t test

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 0.93523, num df = 11, denom df = 11, p-value = 0.9136
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2692306 3.2486944
## sample estimates:
## ratio of variances 
##          0.9352262

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -2.1066, df = 21.975, p-value = 0.0468
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.80197969 -0.00625586
## sample estimates:
## mean of x mean of y 
##  2.263029  2.667146

Run ANCOVA

res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()

## linear fit of dependent variable with covariates

## `geom_smooth()` using formula 'y ~ x'

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6998 -0.9217 -0.2470  0.9448  3.7247 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  29.1574     1.6797  17.358 2.52e-14 ***
## df[[2]]      -0.7788     0.6682  -1.166    0.256    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.614 on 22 degrees of freedom
## Multiple R-squared:  0.05816,    Adjusted R-squared:  0.01535 
## F-statistic: 1.359 on 1 and 22 DF,  p-value: 0.2563
## 
## y = 29 + -0.78x, r^2 = 0.0582

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   12.0000000    0.0000000    0.0000000    1.8745230    3.8776946    2.0031716 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   32.0057573    2.5865853    2.6671464    0.1378974    0.3035102    0.2281884 
##      std.dev     coef.var 
##    0.4776907    0.1791018 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   12.0000000    0.0000000    0.0000000    1.5858365    2.8022937    1.2164572 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   27.1563440    2.3163711    2.2630287    0.1333566    0.2935159    0.2134077 
##      std.dev     coef.var 
##    0.4619608    0.2041338 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.5967 0.4481
##       22               
## _____________________
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  0.980  0.9799   4.438 0.0468 *
## Residuals   22  4.858  0.2208                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.001834 0.0018341   5.682 0.0262 *
## Residuals   22 0.007101 0.0003228                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  0.980  0.9799   4.305 0.0505 .
## df[[3]]      1  0.078  0.0778   0.342 0.5649  
## Residuals   21  4.780  0.2276                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##             Sum Sq Df F value  Pr(>F)  
## (Intercept) 1.0262  1  4.5087 0.04577 *
## df[[1]]     0.7182  1  3.1554 0.09017 .
## df[[3]]     0.0778  1  0.3420 0.56492  
## Residuals   4.7797 21                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

RQ

Note that some datapoints were removed due to measurement error. Data was filtered for values < 1 and > 0.5.

Hide

Select a tab above to view

Total cycle

total.RQ <- total %>% select(contains("RQ")) %>% filter_all(all_vars(. < 1 & . > 0.5))
total.RQ.WT <- total.RQ %>% select_WT()
total.RQ.KO <- total.RQ %>% select_KO()

rcf.WT <- lapply(total.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.0005))
rcf.KO <- lapply(total.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.0005))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in complete cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 0.815989968208129"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 0.820796165978486"

Run F & t test

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.7748, num df = 11, denom df = 11, p-value = 0.3555
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.510926 6.165133
## sample estimates:
## ratio of variances 
##           1.774803

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -0.72408, df = 20.409, p-value = 0.4772
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02871452  0.01390229
## sample estimates:
## mean of x mean of y 
## 0.8160276 0.8234337

Run ANCOVA

res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()

## linear fit of dependent variable with covariates

## `geom_smooth()` using formula 'y ~ x'

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6588 -0.9693  0.2279  0.8568  4.1602 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    17.90      11.30   1.584    0.127
## df[[2]]        11.39      13.78   0.827    0.417
## 
## Residual standard error: 1.638 on 22 degrees of freedom
## Multiple R-squared:  0.03015,    Adjusted R-squared:  -0.01394 
## F-statistic: 0.6838 on 1 and 22 DF,  p-value: 0.4172
## 
## y = 18 + 11x, r^2 = 0.0301

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 1.200000e+01 0.000000e+00 0.000000e+00 7.946162e-01 8.641973e-01 6.958113e-02 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 9.881205e+00 8.237821e-01 8.234337e-01 6.140242e-03 1.351458e-02 4.524308e-04 
##      std.dev     coef.var 
## 2.127042e-02 2.583137e-02 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 1.200000e+01 0.000000e+00 0.000000e+00 7.683244e-01 8.531133e-01 8.478887e-02 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 9.792331e+00 8.192496e-01 8.160276e-01 8.180137e-03 1.800436e-02 8.029757e-04 
##      std.dev     coef.var 
## 2.833683e-02 3.472533e-02 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.7314 0.4016
##       22               
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000329 0.0003291   0.524  0.477
## Residuals   22 0.013809 0.0006277               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 1.087e-05 1.087e-05   3.464 0.0761 .
## Residuals   22 6.902e-05 3.137e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000329 0.0003291   0.531  0.474
## df[[3]]      1 0.000782 0.0007823   1.261  0.274
## Residuals   21 0.013027 0.0006203               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.039605  1 63.8444 8.398e-08 ***
## df[[1]]     0.000685  1  1.1046    0.3052    
## df[[3]]     0.000782  1  1.2611    0.2741    
## Residuals   0.013027 21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

Light cycle

light.RQ <- light %>% select(contains("RQ")) %>% filter_all(all_vars(. < 1 & . > 0.5))
light.RQ.WT <- light.RQ %>% select_WT()
light.RQ.KO <- light.RQ %>% select_KO()

rcf.WT <- lapply(light.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.0005))
rcf.KO <- lapply(light.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.0005))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in light cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 0.772440841371076"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 0.780156358993069"

Run F & t test

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 3.2866, num df = 11, denom df = 11, p-value = 0.06047
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.9461303 11.4165622
## sample estimates:
## ratio of variances 
##           3.286572

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -0.32483, df = 17.127, p-value = 0.7492
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.03087977  0.02263589
## sample estimates:
## mean of x mean of y 
## 0.7740705 0.7781924

Run ANCOVA

res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()

## linear fit of dependent variable with covariates

## `geom_smooth()` using formula 'y ~ x'

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2970 -1.1931  0.2283  1.0840  3.7131 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   11.378      8.168   1.393   0.1775  
## df[[2]]       20.434     10.516   1.943   0.0649 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.537 on 22 degrees of freedom
## Multiple R-squared:  0.1465, Adjusted R-squared:  0.1077 
## F-statistic: 3.776 on 1 and 22 DF,  p-value: 0.0649
## 
## y = 11 + 20x, r^2 = 0.146

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 12.000000000  0.000000000  0.000000000  0.734585496  0.805507815  0.070922319 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  9.338309179  0.777979674  0.778192432  0.006129090  0.013490037  0.000450789 
##      std.dev     coef.var 
##  0.021231791  0.027283472 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   0.71307952   0.85084419   0.13776467 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   9.28884594   0.77121792   0.77407049   0.01111137   0.02445596   0.00148155 
##      std.dev     coef.var 
##   0.03849091   0.04972533 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  2.4256 0.1336
##       22               
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000102 0.0001019   0.106  0.748
## Residuals   22 0.021256 0.0009662               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 8.420e-06 8.418e-06   3.459 0.0763 .
## Residuals   22 5.354e-05 2.434e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq  Mean Sq F value Pr(>F)  
## df[[1]]      1 0.000102 0.000102   0.123 0.7289  
## df[[3]]      1 0.003907 0.003907   4.730 0.0412 *
## Residuals   21 0.017348 0.000826                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##                Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.0231390  1 28.0095 3.015e-05 ***
## df[[1]]     0.0008806  1  1.0659   0.31361    
## df[[3]]     0.0039074  1  4.7299   0.04122 *  
## Residuals   0.0173483 21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

## [1] "Varies significantly with Covar, p = 0.0412152943838807"

Dark cycle

dark.RQ <- dark %>% select(contains("RQ")) %>% filter_all(all_vars(. < 1 & . > 0.5))
dark.RQ.WT <- dark.RQ %>% select_WT()
dark.RQ.KO <- dark.RQ %>% select_KO()

rcf.WT <- lapply(dark.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.0005))
rcf.KO <- lapply(dark.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.0005))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = "RQ", ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in dark cycle", cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 0.844890762938068"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 0.850790371295462"

Run F & t test

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.1212, num df = 11, denom df = 11, p-value = 0.8529
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3227626 3.8946419
## sample estimates:
## ratio of variances 
##            1.12118

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -0.8451, df = 21.928, p-value = 0.4072
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.03766279  0.01585748
## sample estimates:
## mean of x mean of y 
## 0.8454874 0.8563901

Run ANCOVA

res <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar) %>% ANCOVA()

## linear fit of dependent variable with covariates

## `geom_smooth()` using formula 'y ~ x'

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5600 -0.8905  0.0512  0.7171  4.2693 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   26.302      9.403   2.797   0.0105 *
## df[[2]]        1.100     11.042   0.100   0.9216  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.663 on 22 degrees of freedom
## Multiple R-squared:  0.0004505,  Adjusted R-squared:  -0.04498 
## F-statistic: 0.009915 on 1 and 22 DF,  p-value: 0.9216
## 
## y = 26 + 1.1x, r^2 = 0.00045

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 12.000000000  0.000000000  0.000000000  0.828795893  0.937811179  0.109015286 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 10.276680821  0.847528600  0.856390068  0.008857991  0.019496306  0.000941568 
##      std.dev     coef.var 
##  0.030684981  0.035830612 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 12.000000000  0.000000000  0.000000000  0.801022574  0.920923635  0.119901060 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 10.145848961  0.847581860  0.845487413  0.009379353  0.020643818  0.001055667 
##      std.dev     coef.var 
##  0.032491033  0.038428760 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   0.158 0.6949
##       22               
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000713 0.0007132   0.714  0.407
## Residuals   22 0.021970 0.0009986               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 1.357e-05 1.357e-05   3.234 0.0859 .
## Residuals   22 9.231e-05 4.196e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000713 0.0007132   0.687  0.417
## df[[3]]      1 0.000155 0.0001551   0.149  0.703
## Residuals   21 0.021815 0.0010388               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.049979  1 48.1129 7.487e-07 ***
## df[[1]]     0.000858  1  0.8260    0.3737    
## df[[3]]     0.000155  1  0.1493    0.7031    
## Residuals   0.021815 21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

PCA

Libraries

reqpkg <- c("ggfortify","rgl","lfda")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "ggfortify"
## [1] '0.4.10'
## [1] "rgl"
## [1] '0.100.54'
## [1] "lfda"
## [1] '1.1.3'

knitr::knit_hooks$set(webgl = hook_webgl)

Functions

pca_mutate <- function(data, channel) {
  df <- data %>% select(contains(channel))
  WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
  KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
  df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
  
  return(df)
}

plot_pca <- function(pca, data, title="") {
  autoplot(pca, data = na.omit(data), colour = "id", loadings = TRUE, loadings.colour = "black", loadings.label = TRUE, loadings.label.colour = "black", loadings.label.size = 5) + ggtitle(title)
}

Light

Small filtering step: during the day, there were points where the gas analyzers failed, leaving VO2 and RQ values unbelievably low. The rows contianing data for those points were simply discarded.

light <- light %>% slice(-320:-325)

VO2_light <- pca_mutate(light, "VO2")
RQ_light <- light %>% pca_mutate("RQ")
Food_light <- pca_mutate(light, "Food")
Water_light <- pca_mutate(light, "Water")
BM_light <- pca_mutate(light, "BodyMass")
PedMeters_light <- pca_mutate(light, "PedMeters_R")
WheelMeters_light <- pca_mutate(light, "WheelMeters")

df_light <- RQ_light %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_light$value) %>% add_column(FoodIn=Food_light$value) %>% add_column(WaterIn=Water_light$value) %>% add_column(BodyMass=BM_light$value) %>% add_column(PedMeters=PedMeters_light$value) %>% add_column(WheelMeters=WheelMeters_light$value)

PCA of most variables

df_light %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title = "PCA of most variables, light cycle")

PCA without BodyMass

df_light %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA of most variables besides BM, light cycle")

VO2 vs RQ

df_light %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA VO2 and RQ, light cycle")

Food vs Water

df_light %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA Food and Water, light cycle")

Including extra variables (AllMeters, ZBreaks)

AllMeters_light <- pca_mutate(light, "AllMeters_R")
ZB_light <- pca_mutate(light, "ZBreak")
df_light %<>% add_column(AllMeters=AllMeters_light$value) %>% add_column(ZBreak=ZB_light$value)

df_light %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA of all variables, light cycle")

Dark

VO2_dark <- pca_mutate(dark, "VO2")
RQ_dark <- dark %>% pca_mutate("RQ")
Food_dark <- pca_mutate(dark, "Food")
Water_dark <- pca_mutate(dark, "Water")
BM_dark <- pca_mutate(dark, "BodyMass")
PedMeters_dark <- pca_mutate(dark, "PedMeters_R")
WheelMeters_dark <- pca_mutate(dark, "WheelMeters")

df_dark <- RQ_dark %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_dark$value) %>% add_column(FoodIn=Food_dark$value) %>% add_column(WaterIn=Water_dark$value) %>% add_column(BodyMass=BM_dark$value) %>% add_column(PedMeters=PedMeters_dark$value) %>% add_column(WheelMeters=WheelMeters_dark$value)

PCA of most variables

df_dark %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, "PCA of most variables, dark cycle")

PCA without BodyMass

df_dark %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA of most variables besides BM, dark cycle")

VO2 vs RQ

df_dark %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA VO2 and RQ, dark cycle")

Food vs Water

df_dark %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA Food and Water, dark cycle")

Including extra variables (AllMeters, ZBreaks)

AllMeters_dark <- pca_mutate(dark, "AllMeters_R")
ZB_dark <- pca_mutate(dark, "ZBreak")
df_dark %<>% add_column(AllMeters=AllMeters_dark$value) %>% add_column(ZBreak=ZB_dark$value)

df_dark %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA of all variables, dark cycle")

Local Fisher Discriminant Analysis (LDFA)

I attempted to model and plot ldfa. LDFA should give similar results to PCA, with the advantage of being sensitive to multimodality of data, i.e. multiple effects driving an observation, as well as sensitive to distinctiveness in the data that might be lost in reducing dimensions in PCA. In other words, a learning model like LDFA might account for more relationships than a principal component.

More information:

Tang Y, Cohen J, Local Fisher Discriminant Analysis on Beer Style Clustering.
Yang Y, Li W, lfda: An R Package for Local Fisher Discriminant Analysis and Visualization. The R Journal Vol. XX/YY, AAAA 20ZZ

Light

model <- lfda(df_light[-1], df_light[, 1], r=3, metric = "plain")
autoplot(model, data=df_light, frame = T, colour = "id")

plot(model, labels = df_light[, 1])

You must enable Javascript to view this page properly.

Dark

model <- lfda(df_dark[-1], df_dark[, 1], r=3, metric = "plain")
autoplot(model, data=df_dark, frame = T, colour = "id")

plot(model, labels = df_dark[, 1])

You must enable Javascript to view this page properly.

Session info

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lfda_1.1.3       rgl_0.100.54     ggfortify_0.4.10 WRS2_1.0-0      
##  [5] pastecs_1.3.21   multcomp_1.4-12  TH.data_1.0-10   MASS_7.3-51.5   
##  [9] survival_3.1-11  mvtnorm_1.1-0    effects_4.1-4    compute.es_0.2-4
## [13] car_3.0-7        carData_3.0-3    nplr_0.1-7       purrr_0.3.3     
## [17] scales_1.1.0     tibble_2.1.3     docstring_1.0.0  reshape2_1.4.3  
## [21] DT_0.13          ggpubr_0.2.5     ggplot2_3.3.0    dplyr_0.8.5     
## [25] magrittr_1.5    
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.4             colorspace_1.4-1        ggsignif_0.6.0         
##  [4] rio_0.5.16              htmlTable_1.13.3        markdown_1.1           
##  [7] base64enc_0.1-3         mc2d_0.1-18             rstudioapi_0.11        
## [10] roxygen2_7.1.0          farver_2.0.3            RSpectra_0.16-0        
## [13] xml2_1.2.5              codetools_0.2-16        splines_3.6.3          
## [16] knitr_1.28              Formula_1.2-3           jsonlite_1.6.1         
## [19] nloptr_1.2.2.1          cluster_2.1.0           png_0.1-7              
## [22] shiny_1.4.0.2           compiler_3.6.3          backports_1.1.5        
## [25] fastmap_1.0.1           assertthat_0.2.1        Matrix_1.2-18          
## [28] survey_3.37             later_1.0.0             acepack_1.4.1          
## [31] htmltools_0.4.0         tools_3.6.3             gtable_0.3.0           
## [34] glue_1.3.2              Rcpp_1.0.4.6            cellranger_1.1.0       
## [37] vctrs_0.2.4             nlme_3.1-144            crosstalk_1.1.0.1      
## [40] xfun_0.12               stringr_1.4.0           openxlsx_4.1.4         
## [43] lme4_1.1-21             miniUI_0.1.1.1          mime_0.9               
## [46] lifecycle_0.2.0         klippy_0.0.0.9500       zoo_1.8-7              
## [49] hms_0.5.3               promises_1.1.0          sandwich_2.5-1         
## [52] RColorBrewer_1.1-2      yaml_2.2.1              curl_4.3               
## [55] gridExtra_2.3           rpart_4.1-15            reshape_0.8.8          
## [58] latticeExtra_0.6-29     stringi_1.4.6           checkmate_2.0.0        
## [61] manipulateWidget_0.10.1 boot_1.3-24             zip_2.0.4              
## [64] rlang_0.4.5             pkgconfig_2.0.3         evaluate_0.14          
## [67] lattice_0.20-40         htmlwidgets_1.5.1       labeling_0.3           
## [70] cowplot_1.0.0           tidyselect_1.0.0        plyr_1.8.6             
## [73] R6_2.4.1                Hmisc_4.4-0             DBI_1.1.0              
## [76] pillar_1.4.3            haven_2.2.0             foreign_0.8-75         
## [79] withr_2.1.2             mgcv_1.8-31             abind_1.4-5            
## [82] nnet_7.3-13             crayon_1.3.4            rARPACK_0.11-0         
## [85] rmarkdown_2.1           jpeg_0.1-8.1            grid_3.6.3             
## [88] readxl_1.3.1            data.table_1.12.8       forcats_0.5.0          
## [91] webshot_0.5.2           digest_0.6.25           xtable_1.8-4           
## [94] tidyr_1.0.2             httpuv_1.5.2            munsell_0.5.0          
## [97] mitools_2.4

By Sumeed Yoyo Manzoor

smanzoor@uchicago.edu

RCF and PCA of WT vs L-Bmal1-KO RMR wheel data from metabolic cages

Sumeed ManzoorChang Lab at the University of Chicago

July 7, 2020, updated Jul 07, 2020

Dataset

Create data

Import data

Datatables

Light

Hide

VO2

RQ

BodyMass

Food

Water

PedMeters

WheelMeters

Dark

Hide

VO2

RQ

BodyMass

Food

Water

PedMeters

WheelMeters

Graphs

Hide

VO2

RQ

BodyMass

Food

Water

PedMeters

PedSpeed

AllMeters (R)

AllMeters (M)

WheelMeters

WheelSpeed

ZBreaks

H20

VOC

H2

RCF

Libraries

Functions

VO2

Hide

Total cycle

Light cycle

Dark cycle

RQ

Hide

Total cycle

Light cycle

Dark cycle

PCA

Libraries

Functions

Light

Dark

Local Fisher Discriminant Analysis (LDFA)

Light

Dark

Session info

Sumeed Manzoor
Chang Lab at the University of Chicago